function [ar_coeffs,NoiseVariance] = ar_yule(x,p)

[R, lags] = xcorr(x, p, 'biased'); % 计算自相关函数，lags 是滞后值
R = R(lags >= 0); % 只保留非负滞后的自相关函数值
R1 = toeplitz(R(1:p)'); % 构造 Toeplitz 矩阵
r = R(2:1+p); % 构造自相关向量

% 求解 Yule-Walker 方程，得到 AR 系数
ar_coeffs = -R1 \ r; % a_1,a_2,... a_p

% 计算预测误差功率
NoiseVariance = R1(1, 1) + ar_coeffs' * r;
ar_coeffs = [1,transpose(ar_coeffs)]; % same with aryule()
end